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1. Introduction 

Understanding of tlie electronic structure of extended systems with a local perturbation, 
e.g. point defects in the crystal bulk [T] or adsorption of molecules at crystal surfaces 
j2j is of fundamental importance in solid state physics and chemistry. One way of 
calculating the electronic structure of the mentioned systems is based on the usage of 
so-called cluster methods in which a finite fragment of an extended system (a quantum 
cluster) is considered in detail while the rest of the system is treated at a lower level 
of theory [lIlElinilZllHlinilinillllliaill- The main problem of any existing cluster 
based scheme is in choosing an appropriate termination of the cluster. Usually, the 
quantum cluster is surrounded by point charges 0, pseudoatoms (see, e.g. [ISl), link 
atoms [HI El EI or pseudopotentials jTH [13 [12] • In more sophisticated methods the 
environment region is described by an electronic wavefunction which could be either 
frozen ^21 QB] or recalculated self-consistently with that of the quantum cluster region 
[T7| [HI [IHl [ini [2II| (a general theory of cluster embedding which comprises most of the 
existing cluster schemes is considered in jHllll)- 

A rather general cluster method based on overlapping (not orthogonal) localised 
orbitals is presently being developed in our laboratory. Our method which is similar in 
spirit to some one-electron methods [ini [201 [IHl is based on a construction of strongly 
localised orbitals which are designed to represent the true electronic density of the 
entire system via a combination of elementary densities associated in simple cases with 
atoms, ions and/or bonds; these are called regions |2I]. Our intention is to create a 
rather general technique which can be valid for systems of different chemical character, 
ranging from purely ionic to strongly covalent (excluding metals). Therefore, the proper 
choice of the localisation technique as well as a general method of calculating electron 
density out of strongly localised non-orthogonal orbitals localised within corresponding 
regions are crucial for our method to work for a wide range of systems. 

The issue of calculating orbitals localised in appropriate regions for extreme cases 
of strongly ionic and covalent crystals has been considered separately ^Tj. It is the main 
objective of this paper to discuss methods of calculating the electron density of periodic 
systems described via localised non-orthogonal orbitals. 

It should be mentioned that the literature on this topic is quite scarce which is 
probably explained by the lack of interest (until recently) to non-orthogonal (non- 
canonical) molecular orbitals: in most techniques used in the solid state community 
orthogonal Bloch functions are employed in practical calculations. There are only a few 
exceptions (see e.g. |22j). If a set of non-orthogonal orbitals is used, the expression for 
the electron density is much more complicated since it contains an inverse of an infinite 
overlap matrix constructed out of the non-orthogonal orbitals of the whole system under 
consideration j2S]- 

As far as we are aware, there have only been two methods developed which enable 
calculation of the electron density of a periodic system from non-orthogonal orbitals. 
Both methods are based on a series expansion of the density: while the first method 
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Figure 1. A schematic of a crystal division into overlapping regions. Each atom 
(shown by small black circles) provides a set of atomic orbitals centred on it. Only two 
neighbouring regions A (solid line) and B (dashed line) are shown which physically 
represent bonds between atoms 1-2 and 2-3, respectively. All atomic orbitals centred 
on atoms within each region contribute to the localised orbitals associated with this 
region. It is thus seen that different regions may have common atomic orbitals if their 
borders overlap. In particular, atomic orbitals of atoms 1, 2 and 3 belong to both 
indicated regions. 

|22] relies on the so-called cluster expansion of the density, the second one |211 ESI 121] 

is based on the power expansion of the inverse overlap matrix. In this paper we analyse 
only the second of the methods in detail since the first one is very complicated and 
much more difficult to implement. In section 2 we reexamine the second method from 
the point of view of the correct density normalisation. Then, in section 3 we suggest an 
alternative technique which does not require any series expansion. Both methods are 
compared in section 4 using a very simple model system. The paper is finished with a 
short discussion and conclusions in section 5. 

2. Electron density of a periodic system 

Let Capital letters A, B, etc. be used to indicate regions, while the corresponding small 
letters a, b, etc. - localised orbitals associated with them, i.e. a & A, b & B, etc., see 
Fig. ^ Each region may have several localised orbitals. We assume that the orbitals are 
real. They are expanded over atomic orbitals centred only on atoms which are inside the 
region border. Two localised orbitals belonging to different regions are not orthogonal 
either because they have common atomic orbitals or, if they do not, then due to their 
exponential tails. 

Each region A is prescribed with an even number of electrons. Thus, there 
is a finite number = Na/2 of double occupied orbitals associated with the given 
region A. The localised orbitals v'Aa(r) belonging to the same region are assumed 
to be orthonormal; orbitals belonging to different regions are not orthogonal, i.e the 
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corresponding overlap integral SAa,Bb = {y^Aa Ifsb) is assumed to be not zero in general. 
Note that there might be several regions within every primitive unit cell in the crystal. 
Localised orbitals belonging to physically equivalent regions in different unit cells are 
obtained by appropriate lattice translations, i.e. (fsai'''') = V^Aali" — L), where physically 
equivalent regions B and A are separated by the translation vector L. 

In general, the spinless electron density takes on the following general form 

p(r,r') = 2^^v9Aa(r)(S"^)Aa,B6m(r') 



(1) 



Aa Bb 



which contains the inverse of the overlap matrix, S =||5'Aa,_B6|| • The overlap matrix can 
also be written as a set of finite matrix blocks SAB=||S'Aa,Bfe|| associated with every pair 
of regions. Note that for an infinite crystal the matrix S has an infinite size. As usual, 
the factor of two is due to the fact that each orbital is occupied by two electrons with 
opposite spins. 

In both summations above localised orbitals from all unit cells are taken into 
account. To stress the periodic symmetry of the crystal, it is useful to rewrite the 
density in a slightly different form: 

p(r,r')=Ep(r-L,r'-L) (2) 

L 

where we introduced a periodic image of the density 



■'density image" for short): 



P(r,r') = 2^ 'J2vAa{r){S ^) Aa,Bb^ Bb{r' 

Aa Bb 



(3) 



where in the first sum (indicated by a prime) the summation is run only over localised 
orbitals within the single primitive cell associated with the zero lattice translation; the 
other summation runs over all orbitals in the whole infinite system. Note that the 
density image is normalised on the number of electrons in the unit cell only: 

J p{r,r)dr = Y.'NA (4) 



2.1. Method based on the expansion of the inverse of the overlap matrix 



Following the original prescription by Lowdin 



we present the overlap matrix as 



S = 1 + A, where the matrix A =\\A.Aa,Bb \\ is the same as the original overlap matrix 
except for its elements when A = B which are all equal to zero, A.Aa,Aa' = 0. Then, one 
can formally write a matrix expansion: 

S-^= (1 + A)-i = 1- A + A^- A=^ + ... (5) 

One can show (using diagonalisation of the matrix S or its expansion over the 
eigenstates) that the expansion (0) can only be used if absolute values of all eigenvalues 
of the matrix A are less than unity. 

Using expansion of Eq. (jSJ, one obtains the following expansion for the image 
density 0: 



P(r,r') = Y.P 



E(-ir 



re=0 



n=0 



2j2'Y.VAa{r){A'')Aa,Bb^Bb{r') 



Aa Bb 



(6) 
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Note that a general n-th order term (for n > 2) contains additional n — 1 summations 
over all regions due to matrix multiplications in A". 

In principle, formulae (j2I) and © allow an approximate calculation of the electron 
density by terminating the infinite expansion. Care should be taken, however, in doing 
so in order to preserve the correct normalisation of the density. 

The zero order contribution, 

pW(r,r') = 2^ VAa(r)¥.^.(r') (7) 

Aa 

does not contain overlap integrals at all and is normalised to the total number of electrons 
in the unit cell. Therefore, if any higher order terms are kept in the terminated expansion 
(ini), they should be integrated (normalised) to zero. Consider this point in more detail. 
The first order contribution to the image density, 

p«(r,r') = -2^ 'J2vAa{r)AAa,BbVBb{r') (8) 

Aa Bb 

contains different regions A ^ B and thus its contribution to the charge (or 
normalisation) becomes: 

AiV« = / p«(r,r)dr = -2^ 'Y^^BbAa^Aa^Bb = -2^ 'TiA (A^) (9) 

Aa Bb A 

where the trace Tia{- ■ •) here is calculated with respect to the localised orbitals belonging 
to region A only. We see that the first order term has a finite nonzero charge (in fact, 
it is negative). 

Any higher order contributions in Eq. © for n > 2 contains additional summations 
over regions so that equal regions A = B in the double summation there are also possible. 
Therefore, every such contribution, p*^"^(r,r'), will be split into two terms: a diagonal 
term, 

pfV,r ) = 2(-l)" Y: 'VAa{r){\-)AaAa'VAa'{r') (10) 
Aa,a' 

in which A = B, and a non-diagonal term, 

pS(r,rO=2(-l)"E' E ^Aa{r){A-)A^,sb^Bb{r') (11) 

Aa B{j^A),b 

associated with A ^ B in Eq. (jHl). Correspondingly, we obtain the following 
contributions to the charge: 

AiVi") = 2(-l)" 'iAlAa,Aa'SAa,Aa' = 2(-l)"E 'Tr^ (A") (12) 

Aaa' A 

^NirJ = 2i-irY' E (A")Aa,B.ABMa = 2(-l)"E'TrA(A"+i)(13) 

Aa B(^A),b A 

Thus, we see that in any order n > 2 we have = —ANjp~^^\ This means 

that the non-diagonal contribution to the density (jllll is compensated exactly by the 
diagonal one (fTUj) of the next order. For instance, the non-zero charge (jH)) is to be 
exactly eliminated by a charge due to the diagonal second order density; in turn, a 
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nonzero charge due to non-diagonal second order density is compensated exactly by the 
diagonal third order density contribution, and so on. 

This result is very useful since it allows one to balance properly a terminated 
expansion for the image density so that it would correspond (in any order!) to the 
correct total charge. To do this, the final expression for the density of any n-th order 
should also include the diagonal {A = B) term from the contribution of the next order. 
We stress that this fact was ignored in the previous applications of this method P^I26j . 
We obtain, that the correct n— th order expansion for the image density in the notations 
of Eqs. (fTUj). (fTTj) should have the form: 

n 

p(r, r') ^ pN(r, r') ^ [pJV, r') + pS(r, r')] + pf+'V, r') (14) 

1=0 

By employing this termination of the series, the normalisation condition Q is satisfied 
exactly. 

Thus, in order to calculate the density up to the th order, one has to calculate the 
matrix elements (^A^^^ of the matrix A for all powers = 1, . . . , n; in addition, one 

also need diagonal A = B elements of A"+^. Then, the contributions from all density 
images corresponding to all lattice translations, Eq. are to be added together to 
get the final electron density. 

The method described here relies on the convergence of the density expansion 
The better localisation of the orbitals <^yia(r); the faster convergence and thus smaller 
number of terms is needed. We shall demonstrate in section IHl that in some cases of not 
very well localised orbitals one has to consider the density expansion up to a very high 
order which makes the calculation extremely time-consuming. Moreover, if the orbitals 
localisation becomes worse than a certain criteria (to be also discussed in section E}, 
then this method fails altogether as the expansion diverges. A general and an extremely 
efficient technique which is not based on a perturbative expansion of any kind and can 
be used for localised orbitals of practically any degree of localisation is suggested in the 
next subsection. 



2.2. Method based on the Fourier transform of localised orbitals 

In Eq. for the electron density, regions A and B are to be chosen from all unit 
cells of the infinite periodic system. It is convenient in this section to identify explicitly 
the lattice vector for every localised orbital in its index. Therefore, in the following we 
shall use letters A, B, etc. only for regions within the zeroth unit cell; in particular, 
the orbital V5yia(r) is assumed to be from the zeroth cell. Localised orbitals from other 
cells are characterised by the combined index (LAa), i.e. v^Lyia(i') = V'Aa(r — L) is the 
a-th localised orbital from region A in the unit cell separated from the zeroth cell by 
the lattice translation L. 

Correspondingly, Eq. (0) is rewritten in the following way: 

p(r,r') = 2^ J2 <^Aa(r-L)(S-i)LAa,MB6m(r'-M) (15) 

LAa MBb 
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where L and M are two lattice translations. A further transformation is possible here 
since the overlap integral Si^Aa,MBb depends in fact only on the difference M — L of 
the lattice translations. This allows expansion of the overlap integral into the Fourier 
integral 

Sl^AaMBb = ^ E 5Aa,B6(k)e-i^(^-^) (16) 

where the summation is performed over points k in the first Brillouin zone (BZ) and 

SAa,Bb(^) = SoAa,-LBbe^^^ (17) 
L 

is the corresponding Fourier image. The direct lattice summation in the last formula 
is easily terminated due to (usually) exponential decay of the overlap integrals between 
localised orbitals. 

Using the Fourier representation of the overlap matrix, one can exactly calculate 
its inverse as follows: 

fS-M =-yfS-i(k)l e-ii^^^-^) (18) 

k 

Note that the matrix S(k) =|| SAa.Bb^^ || has a finite size of the number of localised 
orbitals per unit cell. Therefore, in order to calculate the inverse of the overlap matrix 
in direct space, one has to perform the calculation of S^^(k) for finite size matrices for 
every k point necessary to sample the BZ. Substituting Eq. (fTHj) into Eq. (fT3j) . we arrive 
at the following final expression for the electron density: 

P(r,r') = |E|EE^Aa(r,k) [S-^(k)l^^^^y,^,(r', k)| (19) 



^ k y Aa Bb 



where 



^Aair, k) = 5] ^Aaij " L)e-i'^^ (20) 

L 

is the Fourier expansion of the localised orbital. Due to exponential decay of the localised 
orbitals, the summation over lattice vectors L in the last expression is in fact finite. 

The obtained formula for the density is exact. In particular, it contains the 
periodicity of the lattice built in. It is also extremely convenient for numerical 
implementation. Indeed, what is needed is the calculation of the Fourier images, 
according to Eq. ()20|) . of every localised orbital in the primitive unit cell for every 
k point. The summations in the curly brackets in Eq. (fTTH) are finite (limited to the 
orbitals within the zeroth cell only) and are thus easily performed. The extend to 
which the orbitals (y9yia(r) are localised is reflected by the number of cells to be taken 
into account while performing the lattice summations in Eqs. (jl7p and (j2(Jj) . Even for 
orbitals which are not very well localised, the amount of work needed to perform these 
lattice summations is incomparable with the cost of the first method (section 12.11] which 
requires including more terms in the perturbation expansion if the localisation is not 
good enough. 
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3. Results 

Atomic units are used throughout this section. The apphcation of the two methods 
considered in the previous sections is illustrated here on a simple cubic lattice model 
containing a single region in every unit cell. The lattice constant a will be assumed to 
be equal to 1 a.u. for simplicity. Each region is represented by a single localised orbital 
in a form of a normalised s type Gaussian 

V^LAa(r) - ^L(r) = ^(r - L), ^(r) = (^-j e'"^ . (21) 

By choosing various values for the exponent a, one can vary the degree of localisation 
of the orbitals. Indeed, the size of the spatial extent of the orbital can be measured in 
terms of r^jf = \f^^ — 1.52a^^/^, which corresponds to e~°'^'^ff =0.1. We found this 
approach more convenient in our particular case than the application of the existing 
localisation criteria (see, e.g. j2Zll2H])- 

For this model system it is possible to do some preliminary analytical estimates of 
the convergence of the series ©. We know from section ITT] that the series will converge 
if all eigenvalues of the matrix A = S — 1 are between -1 and 1. It is easy to 
notice that the eigenvalues are in fact given by the Fourier transforms of the matrix 
A which is introduced much in the same way as S(k) in Eq. |T7|). Indeed, because 
Alm = AoM-L, one can write: 



M \ M / 



(22) 



This is nothing but the eigenproblem for the matrix A with A^ being its eigenvalues 
(numbered by vectors k from the BZ) and e^^^ - eigenvectors. Therefore, the 
convergence criteria for the series @ reduces to the inequalities |Ak| < 1 which should 
be valid for any k. Taking into account the overlap only between nearest neighbours, 
we obtain: 

|Ak| = \26{cos{kxa) + cos{kya) + cos(A;2a))| < 65 < 1 

with the overlap between neighbouring orbitals being 6 = e~°"^/^. This results in the 
following criterion for the convergence of the Lowdin expansion (for a =1 a.u.): 

ayal = 21n(6) ^ 3.6 (23) 

Similar analysis which takes into account the next nearest neighbours gives a very similar 
estimate of al ~ 4.05. These estimates correspond to the maximum spatial extent of 
the orbitals ()2H1 of the order of Veff ~0.76 a.u., i.e. there is very small overlap between 
neighbouring orbitals which, we recall, are separated by 1 a.u. in the lattice. 

The other method based on the Fourier transform of the orbitals has also its limits 
which are hidden in the formulae ()17j) and ()20|1 : if a certain cut-off |L| < Vc for the 
direct lattice summation L is assumed in the calculation of SAa,Bb(X) and (/9yia(r, k), 
then there will be some limitations on the allowed degree of localisation of the orbitals. 
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The required criterion can be worked out e.g. by analysing the Fourier transform (|2(jp of 
the orbital at its maximum in the centre of the BZ (i.e. of (p{r = 0, k = 0)) as follows: 

|L|>rc |L|<rc 

Replacing the sums by the corresponding volume integrals, we obtain the following 
criterion: 

xe + - — erfc(x) < - — 

where x = Tcy/a. The inequality above is satisfied if x ^ 2, i.e. a y Assuming 
that Tc is equal to 4-h5 lattice constants, we obtain the necessary condition for the 
exponent of the localised orbitals, 

a > ~ 0.2 (24) 

for which our Fourier transform method should work. The obtained critical value of 
results in the maximum spatial extent of the orbitals of the order of rg// ~3.4 a.u. 
which corresponds to very diffuse orbitals spreading over more than six unit cells. 

Similar criteria is obtained for the ovelrap integrals as well. Thus, the method we 
suggest should have a much wider range of applicability than the Lowdin method as far 
as the degree of localisation of the non-orthogonal orbitals is concerned since ^ ctj. 
This conclusion is also supported by our numerical calculations which we now describe. 

Numerical calculations of the necessary powers of the A matrix needed for the 
Lowdin method were done in the following way. Since the density is calculated in the 
same point r =r' in Eqs. ()10j). and ()14|) . the regions A and B in these equations 
are either the same or not far away from each other. Therefore, to calculate ( A")^^ 
one can simply choose a sufficiently big finite cluster of atoms (in fact, the cluster radius 
should be at least of the order of |r*, where r* is the decay length of the overlap integral) 
with regions A and B somewhere in its centre and then calculate the complete ovelrap 
matrix for it, A. Then, by performing the necessary n — 1 matrix multiplications, one 
can calculate (A") 4„ ri, as ( A" ) 

V ^Aa,Bb V / Aa,Bb 

When using the Fourier transform method, we employed the Monkhorst-Pack (MP) 
method [23 for the k point sampling and the same cut-off distance for the direct lattice 
summations in Eqs. (|T7jl and (j20|) as in the previous method. In all our calculations we 
used the 4x4x4 MP set which was found to be sufficient in all cases. 

Results of our calculations for a large value of the exponent (a ^ ^ 0:2) are 
shown in FigEl This case corresponds to strongly localised orbitals as is the case in ionic 
systems such as MgO and NaCl. Overlap between orbitals is negligible and even zero 
order approximation Lowdin method, Eq. ((Zj), was found sufficient to give the correct 
density. Density curves for both methods are indistinguishable from each other. 

The calculated densities in the intermediate case [a ~ aj) are shown in Fig. El 
This value of a may correspond to ion-covalent and covalent systems. One can see that 
high order approximations (up to n =8) of the Lowdin method, Eq. (HD), are needed 
here to converge the density and thus the calculation is quite time consuming. 
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Figure 2. The exact electron density p{r,r), Eq. I|19|l . and that based on the zero 
order approximation p'^"^ (r, r) , Eq. {Tj), both calculated along the (100) direction using 
a =10 a.u. Note that the densities are nearly zero between the localisation centres 
shown by grey circles. 




Figure 3. The exact electron density p(r, r) , Eq. I|19|) . and the several approximations 
to it using Eq. (|14|l with n —0, 1, 3 and 8, all calculated along the (100) direction 
using a =4 a.u. Note that the density is small (but nonzero) between the localisation 
centres. 
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Figure 4. The electron densities for a =2 a.u. Other notations are the same as in 
Fig. 121 Note that the correct density (dots) is rather large between the localisation 
centres. 

Finally, we show in Fig. EJtlie densities calculated using both methods for orbitals 
which are least localised when a* > a > The density, obtained using the Fourier 
transform method, Eq. (fT^ . is spread almost uniformly in the crystal volume and thus 
may correspond to a metallic band. At the same time, the Lowdin expansion method, 
Eq. dm), does not converge at all and the density is clearly diverges. 

One can expect that the latter situation can happen only for metallic systems. 
Interestingly, our calculations (not reported here) for such a realistic covalent system as 
a crystalline Si show that the Lowdin approach also fails in some cases when the orbitals 
are not sufficiently localised. Note that various degree of localisation of the orbitals can 
be obtained using different localisation techniques and different choice of regions, see 
fIT\ for more details. 

4. Conclusions 

In summary, we have considered two numerical methods which allow calculation of 
the electron density of a 3D periodic system constructed via a set of non-orthogonal 
localised orbitals. The first, so-called Lowdin, method based on the power expansion 
of the inverse of the overlap matrix has been found to be efficient only for strongly 
localised orbitals. For an intermediate degree of orbitals localisation this method has 
been found to be quite computationally demanding since many terms in the series are 
to be retained. However, if orbitals are not sufficiently localised (the exact criterion 
has also been suggested), the method fails altogether and the power expansion has been 
shown to be divergent. 



Calculation of electron density of periodic systems. 



12 



Then, we have suggested another method based on the Fourier transform of the 
locahsed orbitals which involves calculations of inverse of only finite matrices and a k 
point summation over the Brillouin zone. This method is computationally much less 
demanding and does not have any convergency problems. Using a simple model for the 
crystal electron density represented via a set of Gaussian s type orbitals in a simple 
cubic lattice (one orbital per unit cell), we have shown that our method works equally 
well within a rather wide range of orbitals having different localisation, whereas the first 
method fails for a relatively weakly localised orbitals. The application of the Fourier 
transform method to realistic systems such as MgO and Si perfect crystals is published 
elsewhere fIT\ . 
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